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Abstract 



The classical numerical methods play important roles in solving 
wave equation, e.g. finite difference time domain method. However, 
their computational domain are limited to flat space and the time. This 
paper deals with the description of discrete exterior calculus method for 
• numerical simulation of wave equation. The advantage of this method 

is that it can be used to compute equation on the space manifold 
and the time. The analysis of its stable condition and error is also 
| accomplished. 
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X ' 1 Introduction 
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The wave equation is the prototypical example of a hyperbolic partial differ- 
ential equation of waves, such as sound waves, light waves and water waves. 
It arises in helds such as acoustics, electromagnetism, and fluid dynam- 
ics [H[2] . To investigate the predictions of wave equation of such phenomena 
it is often necessary to approximate its solution numerically. A technique 
suitable for providing numerical solutions to the wave propagation problem 
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is the finite difference time domain (FDTD) method. This is normally de- 
fined by looking for an approximate solution on a uniform mesh of points 
and by replacing the derivatives in the differential equation by difference 
quotients at points of mesh. The computational domain of this algorithm is 
limited to flat space and the time [3H7]. 

Discrete exterior calculus (DEC) constitutes a discrete realization of the 
exterior differential forms, and therefore, the right framework in which to 
develop a discretization for differential equations not just on flat space but 
on manifold [8til6|. The differential operators such as gradient, divergence, 
and Laplace operator on manifold can also be naturally discretized using 
DEC. The numerical solution of wave equation on space manifold and the 
time by the methods of DEC is obtained in this paper. For this equation, an 
explicit scheme is derived. The analysis of this scheme's stability shows that 
the numerical solution becomes unstable unless the time step is restricted. 

2 DEC method for wave equation 
Wave equation 

The wave equation is the prototypical example of a hyperbolic partial dif- 
ferential equation. In its simplest form, the wave equation refers to a scalar 
function u that satisfies: 



where A is the Laplace operator and c is the propagation speed of the 
wave. More realistic differential equations for waves allow the speed of wave 
propagation to vary with the frequency of the wave, a phenomenon known 
as dispersion. In this case, c is replaced by the phase velocity: 



Another common correction in realistic systems is that the speed is depend 
on the amplitude of the wave, leading to a nonlinear wave equation: 



DEC scheme for wave equation 

A discrete differential /c-form, k € Z, is the evaluation of the differential k- 
form on all /c-simplices. A dual form is evaluated on the dual cell. Suppose 
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each simplex contains its circumcenter. The circumcentric dual cell D{oq) 
of simplex <jq is 

D(cr ) := (J lnt(c(a )c(ai) ■ ■ ■ c(a T )), 

CToGo-ig-'-ecr,. 

where Oi is all the simplices which contains (To,..., <J%-\, c(o~i) is the cir- 
cumcenter of (jj. In DEC, the basic operators in differential geometry are 
approximated as follows: 

1. Discrete exterior derivative d, this operator is the transpose of the 
incidence matrix of fc-cells on k + 1-cells. 

2. Discrete Hodge Star *, the operator scales the cells by the volumes of 
the corresponding dual and primal cells. 

3. Discrete Laplace operator is *~ 1 d T * +d T * d. 

For some situations, a source having azimuthal symmetry about its axis 
is considered. In this case, the 2D triangular discrete manifold as the space 
is only need to be considered. Now, we show how to derive an explicit DEC 
scheme for Eq.(l) in 2D space manifold and the time. The wave equation 
in 3D space and the time can also be computed by a similar approach. 

Take Fig.l as an example for a part of 2D mesh, in which 0,..., F are 
vertices, 1,...,6 are the circumcenters of triangles, a,...,f are the circumcen- 
ters of edges. Denote hj as the length of line segment and Aijki as the 
area of quadrangle k, I). 
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Fig.l. A part of 2D mesh 

Define 

I12 '■= hf + hf) hz '■= ha + ha, hi '■= he + he, 
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and 

-Pl23456 : = A)l/e + A 02 f a H h A)6de- 

The diffusion term Au at vertice is approximated using discrete Laplace 
operator as follows: 



Au « — ( - u ) + ^(u B - u ) + ^-(u c - u ) 

^123456 V'AO I'BO 'CO 

+ (UD - U ) + (U E ~ U ) + -, [U F - U , 

'DO '£0 'FO 



(2) 



The temporal derivative is approximated by middle time differences as fol- 
lows: 

pp. n i 

-aWW? {t '" + '- 2u " + u "^ (3) 

where At is uniform spacing, and nAt is the coordinate of time. The ap- 
proximation of Eq.(l) generated by substituting the left-hand sides of (2) 
and (3) into (1), thus satisfies 

Right(2)"- 1 = «-2<- 1 + <- 2 ). (4) 

3 Stability, convergence and accuracy 
Stability 

The Courant-Friedrichs-Lewy condition is a necessary condition for stability 
while solving certain partial differential equations numerically. Now, this 
condition is derived for scheme (4). First, this DEC scheme is decomposed 
into temporal and spacial eigenvalue problems. 
The temporal eigenvalue problem: 

d 2 u n . _ 
It can be approximated by difference equation 



(Atf 

Supposing 



A< (5) 



u% +L = t# cos(At) u^' L = cos(-At) 
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and substituting those into Eq.(5), we obtain 

cos(At) + cos(-At) - 2 

therefore 

4 

< A < 0. 



(At)*- " 

This is the stable condition for the temporal eigenvalue problem. 
The spacial eigenvalue problem: 

c 2 Au = Au 

It can be approximated by difference equation (6) based on Fig.l. 



-Pl23456 a ^23 , -.,^34, x . ^45 , 

3 — Au = - — {u A -u ) + - — {u b -uq) + - — (uc-uq) 

■■od i \ , 'lb / \ . L 12 i \ 



+ -, (UD ~ Uq) + 7 (W£ - Uq) + (U F ~ Uq) 

''DO IE0 *F0 

Let Ui = UQCOs(cloi) and substitute into Eq.(6) to obtain 

-^123456 . ^23/ r i \ i \ , ^34 , / , \ n x . ^45 / / ; \ i \ 

5 — A = - — (cos(cZoa) - 1) + -, — (cos(dos) - 1) + — (cos (doc J - 1) 

C lA(j t-B(j 'CO, 

+-^-(cos(d £)) - 1) + t^-(cos(c1 0E ) - 1) + T^-(cos(d F) - 1) 

t_D0 «E0 'F0 

So we have 

2c 2 / ^23_ ^34_ ^45_ _^56_ _^16_ ^12 

-Pl23456 VAO ^50 ^C0 ^D0 ^F0 ^F0 



In order to keep the stability of scheme (4), we need 

4 < 2c 2 / ^ ^ hi_ _|_ _^45_ _|_ _^56_ _|_ _fl6_ _|_ ^ 



(At) 2 -Pl23456 VAO ^B0 ^C0 IdO ho IfO ) 

for all vertices, namely 



cAt < Minoev 



2Pi 



123456 



23 _^34_ _^45_ _^56_ _^16_ |l2_ 
\ ^A0 ^B0 ^C0 ^D0 ^F0 ^F0 
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Convergence 

By the definition of truncation error, the solution u of the Eq.(l) satisfies 
the same relation as scheme (4) except for an additional term O(At) 2 on 
the right hand side. Thus the error Xf = itf — uf is determined from the 
relation 

vn — OV"" 1 yn-2 , (cAt) 2 ( hz , vn -l yn -l\ 

A23456 V'AO 

+ ^(X^ 1 - X™- 1 ) + ^(X^ 1 - X^ 1 )) + O(At) 2 . 

(8) 

Define 

\X n \ = Maxigvl-Xfl- 
From condition (7), we have 

(cAt) 2 / ^23_ _^34_ _^45_ _^56_ _^16_ ^12 

-Pl23456 VAO ^BO ^CO ^DO ^EO ^FO 

Hence the coefficient of Xq _1 in Eq.(8) is nonnegative. It follows that 
\X%\ < 2|A' n - 1 | + \X n ~ 2 \ + 0{At) 2 , 

and hence that 

\X n \ < 2IX™- 1 ! + \X n ~ 2 \ + 0{At) 2 . 
Iterating n, we obtain 

\X n \ < M x \X l \ + M \X°\ + 0(At) 2 , 

where Mi, Mq are finite value defined on n. Since the initial conditions 
ensure X° = and X 1 = 0, we have 

lim \X n \ = 0. 
At->0 

That is to say the numerical solution approaches the exact solution as the 
step size goes to 0, and scheme (4) is convergent. 
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Accuracy 

In scheme (4), the space derivative of is approximated by first order dif- 
ference. Equivalently, u is approximated by linear interpolation functions. 
Consulting the definition about accuracy of finite volume method, we can 
say that scheme (4) has first order spacial accuracy. Scheme (4) has second 
order temporal accuracy, and second order spacial accuracy on rectangular 
grid with uniform spacing. 

4 Algorithm Implementation 

The implementation of DEC scheme (4) for wave equation consists of the 
following steps: 

1. Set the simulation parameters. These are the dimensions of the com- 
putational mesh and the size of the time step, etc.; 

2. Initialize the mesh indexes. 

3. Assign current transmitted signal. 

4. Compute the value of all spatial nodes and temporarily store the result 
in the circular buffer for further computation. 

5. Visualize the currently computed grid of spatial nodes. 

6. Repeat the process from the step 3, until reach the desired total num- 
ber of iterations. 
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The flowchart of the scheme (4) can be seen in Fig. 2. 



Start 



Set propagating media parameters 



Initialize grid indexes 



Assign current transmitted signal 

I 




no 








Store computed grid 




Fig. 2. The flowchart of scheme (4) 

The Fig. 3 and Fig. 4 show the numerical simulation of Gaussian pulse 
propagating on the sphere and rabbit by scheme (4) on C# platform. 
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Fig. 3. The propagation of Gaussian pulse on a rabbit 
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Fig. 4. The propagation of Gaussian pulse on a sphere 

5 Discussion 

The DEC scheme for Laplace operator here can also be used to simulate the 
heat equation, Laplace equation and Poisson equation on manifold. 

Discrete Laplace equation 

The discrete Laplace equation on surface of regular tetrahedron (Fig. (5)) is 
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(9) 



The solution of Eq.(9) is 

ua = ub = uc = u D = C, 

where C is arbitrary constant. Eq.(9) is an imprecise approximation of 
Laplace's equations on a sphere. Obviously, this equation has constant so- 
lution. 



A 




C 

Fig. 5. The surface of regular tetrahedron 
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Discrete Poisson equation 

Consider a discrete Poisson equation on surface of regular tetrahedron. Sup- 
pose the boundary condition is 

u A = H, 

then discrete Poisson equation on Fig. (5) is 



(10) 
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The solution of Eqs.(10) is 



u B = u c = u D = H. 



Discrete heat equation 

The heat equation of temperature u is 

du 



dt 



cAu, 



which can be approximated as 



i — u 



ifi- 1 + cAt Right (2 N 



n-l 



fill 



The Fig. 6 shows the heat diffusion of a constant heat source on the sphere 
simulated by scheme (11). 





Fig. 6. The heat diffusion on a sphere 
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